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We consider the statistics of time-integrated energy fluctuations of a driven bosonic resonator 
(as measured by a QND detector), using the standard Keldysh prescription to define higher mo- 
ments. We find that due to an effective cascading of fluctuations, these statistics are surprisingly 
non-classical: the low-temperature, quantum probability distribution is not equivalent to the high- 
temperature classical distribution evaluated at some effective temperature. Moreover, for a suffi- 
ciently large drive detuning and low temperatures, the Keldysh-ordered quasi-probabihty distribu- 
tion characterizing these fiuctuations fails to be positive-definite; this is similar to the full counting 
statistics of charge in superconducting systems. We argue that this indicates a kind of non-classical 
behaviour akin to that tested by Leggett-Garg inequalities. 

PACS numbers: 



I. INTRODUCTION 

The statistics of photon fluctuations in various setting 
is by now an almost textbook phenomena. Most familiar 
are the statistics that would be measured by a photode- 
tector. Relatively less attention has been paid to photon 
fluctuations in the case where the detection is done in a 
non-demolition manner, meaning that energy quanta are 
measured without destroying them. The study of such 
fluctuations is not just a theoretical curiousity, as quan- 
tum non-demolition (QND) photon detection is experi- 
mentally feasible both optically [T] , as well as in both cav- 
ity QED systems [5H3] and superconducting circuit-QED 
systems .5,, 6]. In the latter systems, one uses disper- 
sive interactions to detect photon number inside a cavity. 
QND detection of phonon number in a mechanical res- 
onator may also soon be possible in optomechanical sys- 
tems [3 [S] , where the energy of a mechanical resonator 
is directly coupled to the frequency of an optical cavity. 

Motivated by developments in optomechanics, we re- 
cently investigated the low-frequency energy fluctua- 
tions of a driven, damped harmonic resonator, focus- 
ing on the possibility of measuring these fluctuations 
non-destructively using an optomechanical cavity [9'. In 
the zero-temperature, quantum limit, the instantaneous 
state of such a resonator is simply a coherent state, yield- 
ing a Poissonian distribution of phonon number. Our 
focus was instead on understanding how the mechanical 
phonon number h fluctuated in time. These fluctuations 
are characterized by a power spectral density S'„ri[w], or 
equivalently by the second central moment of the time- 
integrated phonon number to: 

TO = / dt'n{t'). (1) 

"'0 

As could be anticipated, both S'„„[aj] and {{Srfi)'^) — 
{m^) — {m)'^ have a low- frequency "shot-noise" term pro- 
portional to the average number of phonons induced by 
the drive, n^r] detecting this shot noise contribution 
would be direct evidence for the quantization of the me- 
chanical resonator's energy. As the quantum signature 



here scales as ndr ^ 1, measuring these low- frequency 
energy fluctuations is an easier way of detecting quan- 
tum behaviour than attempting to resolve the instan- 
taneous phonon number and individual quantum jumps. 
Our study also addressed the non-Gaussian nature of the 
driven energy fluctuations by calculating the third cen- 
tral moment {(Srh)^); surprisingly, we found that while 
this quantity is always positive classically, it could be- 
come negative in the low-temperature quantum limit. As 
such, the third moment is far more sensitive to classical- 
quantum differences than the second moment. 

To fully understand the significance of this result, one 
needs to consider the full probability distribution char- 
acterizing the low-frequency fluctuations of n, and com- 
pare its form in the classical and quantum limits. This is 
the objective of this paper: we calculate the distribution 
P{m) in the long-time limit using the standard Keldysh 
operator ordering [TUHH]. Wc find that the anomalous 
negative value of {{5rhY) results from a kind of cascaded 
fiuctuation effect [T31 E] , which can be heuristically at- 
tributed to a correlated fluctuation in the resonator tem- 
perature. We also find that this negative skewness is a 
precusor of something rather dramatic: in the quantum 
limit, the fluctuations of to are most naturally described 
by a quasi-probability distribution P{m) which is not pos- 
itive definite. Such negative counting statistics have been 
encountered before in the study of charge transfer in su- 
perconducting systems [T^ [TS]; their interpretation re- 
quires some care. As we discuss in some detail, they are 
indicative of non-classical temporal correlations, and are 
thus somewhat similar to having violated a Leggett-Garg 
inequality [IB]. Detecting these effects thus represents a 
new way of detecting non-classical behaviour in a driven 
quantum resonator. 

The remainder of this paper is structured as follows. 
In Sec. [nj we introduce our basic model, and present 
our main results for the generating function of P(to); we 
also give a compact review of the Keldysh ordering of 
higher moments for those not familiar with this topic. In 
Sec. |III[ we discuss the form of the distribution in the 
classical limit. Sec. llVl is dedicated to the distribution in 
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the quantum limit, while Sec.[V]is devoted to interpreting 
the negative quasi-probabilities which emerge. Finally, in 
Sec. IVII we discuss issues related to the measurement of 
these effects. An appendix is included which shows how 
the Keldysh operator ordering emerges naturally in the 
proposed experimental realization of Ref. [9], where m 
is measured by using homodyne interferometry to detect 
the frequency shift of an auxiliary cavity. Finally, we note 
that a fermionic analogue of the present problem, the full 
counting statistics of electronic charge fluctuations in a 
chaotic quantum dot, were studied in Ref. [IT] . 

II. MODEL AND CALCULATION 

A. Statement of the problem 

Our damped, driven harmonic resonator is described 
by the Hamiltonian: 

H = Ho + H-y^ hujMc'c - hf (e'""*c + /i.e.) + H^. (2) 

Here, the first term describes the resonator (frequency 
(jJm, number operator h = c^c), describes the damp- 
ing (at a rate 7) and heating of the oscillator by a thermal 
bath, and / is the magnitude of the coherent oscillator 
driving force (frequency lod ^ + S). We take to 
correspond to the standard model of a linear coupling to 
an Ohmic oscillator bath. We also define the dimension- 
less oscillator force susceptibility as: 

^ " 1+4(5/7)2- 

We will be interested in the statistics of the time- 
integrated energy rh (c.f. Eq. ([T])) in the case where 
the oscillator has equilibrated to both the driving force 
and thermal bath long before the initial time t = 0. 
We will also focus exclusively on the long time limit, 
e.g. an integration time t which is long compared to 
1/7. The average and second moment of m are easily 
found by solving the Heisenberg-Langevin equations for 
our system [18l[l9]. In the long time limit, the average 
{rh) ~ [n^-c + nth] where fith denote the thermal num- 
ber of oscillator quanta (determined by the bath tem- 
perature), and fidr — (2//7)^x is the average number of 
quanta due to the driving force. 

For the second central moment, we find in the long 
time limit: 

{{5mf) ^ {{5mY)^r + {{5rnf)t^, (4) 

where 

ux-\2\ 2nth(l + "-th)t i^-s 
{{5m) )th = (5) 

7 

represents a purely thermal contribution whereas 

{{SrnfU. - ^ {n.. + ^) (6) 



represents extra energy fluctuations due to the driving 
force. The last term in Eq.([6| here survives in the limit of 
zero temperature (i.e. nth 0), and is a quantum effect: 
it corresponds to the shot noise fluctuations arising from 
the discreteness of the resonator's energy. 

B. Higher moments and the Keldysh ordering 

Before calculating higher moments and the full distri- 
bution of rh, we must pause to consider the operator- 
ordering ambiguity arising from the non-commutativity 
of n{t) at different times. In calculating the second mo- 
ment, we have naively defined the variance as {fn?) , an 
expression which is naturally symmetrized in terms of n 
products, i.e.: 

(m^) = [ dti [ dt2{h{ti)h{t2)) 
Jo Jo 

= I f dh f dt2{{h{h),h{t2)}). (7) 

^ Jo Jo 

If we define higher moments in the same way (e.g. define 
the jth moment to be {rh^)), they too would be natu- 
rally symmetrized; one might expect that this is then a 
sensible way to proceed. 

Unfortunately, being sensible is not enough to guar- 
antee physical relevance: similar to the standard theory 
of photodetection [50], one must instead model the ac- 
tual detection scheme to properly understand how the 
measured moments correspond to a given operator or- 
dering. As we are interested in non- destructive detection, 
the answer here will not be the normal-ordering prescrip- 
tion used in photodetection. A similar problem arises in 
the measurement of current fluctuations in quantum co- 
herent conductors; the answer emerging from studies of 
this question is the so-called Keldysh operator ordering. 
This ordering appears naturally in a number of idealized 
measurement setups [TDHH]; it can also be given an ele- 
gant motivation using a path-integral formulation of the 
Keldysh field-theoretic technique [5TJ [32] . For the second 
moment, the ordering coincides with the simple definition 
in Eq. Q; for higher moments, the ordering prescriptions 
have no simple intuitive form (see, e.g., Ref.|5] for the ex- 
plicit form for the third moment). 

C. Using an auxiliary qubit to obtain P(m) 

We give here a quick derivation of the Keldysh order- 
ing, and use it to derive P{m) for our system. Following 
Ref. |llj , we consider an idealized method for measuring 
P{m), in which h couples dispersively to the Uz operator 
of an auxiliary two-level system (TLS) with a coupling 
strength fc/2: 

-Hint = -^ncr^. (8) 
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As there are no other terms involving the TLS in the 
Hamiltonian, we see that it simply experiences a mag- 
netic field (X n. If n were just a classical, time-dependent 
field n{t), then during the time interval between and t, 
the TLS would precess an angle d — k J* n{t')dt' = km. 
If now m was a classically stochastic variable described 
by the distribution P{m), then the average of e~*^ over 
this distribution (at a fixed coupling k) would be: 



dmP{m)e 



-ikm 



A[fc]. 



(9) 



Thus, when viewed as a function of fc, the average of the 
precession phase directly yields the moment generating 
function A[fc] of the distribution P{m). 

The above correspondence now provides a means for 
defining P{m) in the quantum case 11 : we simply use 
the fact that the average on the LHS of Eq. ^ corre- 
sponds to p^^(t)/p^^(0), where pi;^{t) is an off-diagonal 
matrix element of the TLS's reduced density matrix. We 
can thus define the moment generating function A[fc] (and 
hence P[m)) in the quantum case via: 



A[fc] 



m(o) 



Tr 



sys 



U{t;k)p,y, (U{t;^k) 



,(10) 



where the time evolution operator U is defined as: 



U(t; k) ^ T exp 



.(11) 



Here, H is given in Eq. 
system (i.e. cavity plus 



(2|, Psys is the initial measured 
Dath) density matrix, and the 
trace is taken over all system degrees of freedom. Further, 
the symbol T denotes time ordering. Eq. ( 10 ) uniquely 



specifies the operating ordering to use for each moment 
of P{m); this is the Keldysh ordering. We stress that 
the same ordering emerges in the analysis of other ideal- 
ized measurement setups [12j : we also show in Appendix 
A that it applies to a realistic setup where n is coupled 
dispersively to a detector cavity whose frequency is mon- 
itored using homodyne detection. 

In our case, the above scheme not only motivates the 
Keldysh ordering, it also gives us a convenient way to 
calculate the generating function A[fc]. The reduced den- 
sity matrix p describing the TLS and the driven resonator 
(i.e. only the resonator's dissipative environment is traced 
out) obeys the following standard master equation: 



Ho,p 



7(nth + l)2?[c]p + 7nth2?[c^]p(12) 



where Hq is defined in Eq. ^ , and where for any operator 

A we define V[A]p ^ ApA^ - (^A^Ap + pA^A^ /2. We 

are using the "quantum optics" version of the master 
equation, which is appropriate for the high-Q limit we 
consider. 

As shown in Ref . jJS] , using standard phase space tech- 



niques, one can solve Eq. (12 1 and thus directly obtain 



A[k]. Ref. [23] used this quantity to study dephasing and 
coherence revivals of the TLS; the emphasis was on un- 
derstanding P'f-iit) for a fixed value of the coupling k. 
In contrast, our focus here is on how p-\i{t) behaves as a 
function of k in the long-time limit, as it is this behaviour 
which will determine P{m) in the long-time limit. 



D. Generating function for P(m) 

Using the procedure described above, and taking the 
long-time limit, the final result for the moment generat- 
ing function A[k] has the simple form (c.f. Eq. (22) in 
Ref. [23]): 



A[fc] - AdrWAthW 



(13) 



where Ath[^] describes a purely thermal (drive- 
independent) contribution, and Adi [fc] describes addi- 
tional fluctuations related to the drive. One finds: 



Adr[fc] 

Ath[fc] 



exp 



= exp 



-ikn^rt 



1 + 4«x (nth + I) (fc/7) - X{k/iy 



,(14a) 



7 7 



(14b) 



In taking the long-time limit, we have simply dropped 
terms in A[k] which decay exponentially in time as 
exp(— 7^/2) or faster. 

Recalling that the jth cumulant of m is give n by 



lnA[fc] 



k=0 



one can easily check that Eqs. ( 14a 



and ( 14b ) yield the same second and third cumulants ob- 



tained from the Heisenberg-Langevin approach. We see 
that the purely thermal fluctuations described by Ath are 
independent of the additional drive-induced fluctuations 
described by A^r; further note that these purely ther- 
mal fluctuations vanish in the limit of zero temperature. 
Thus, in the remainder of the paper we focus on the case 
n-dr ^ "-th, 1, and thus focus attention to the distribution 
P(m) generated by Adr[fc]- 



III. ENERGY FLUCTUATION STATISTICS IN 
THE CLASSICAL LIMIT 

To gain some intuition, it is useful to first consider the 
driven energy counting statistics in the classical, high- 
temperature limit. Formally, one transforms the distri- 
bution described by Eq. ( 14a[ ) to a distribution describing 
the time- integrated energy s = /J dt' E{t') = huiMtn, One 
can then rigorously take the ?i — > limit. Transforming 
back to our original variable to, one finds: 



Adr,ci[fc] = exp 



-ikn^-ct 



1 -f 4ixnth(fc/7) 



(15) 
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In the long time-limit of interest, the corresponding prob- 
ability distribution function can be found within a saddle 
point-approximation, yielding: 



1 



( 



with 



• exp 



o-ci = 2nthX/7- 



, (16) 



(17) 



We see that classically, m is well-approximated as being 
the square of Gaussian random variable with mean 
and standard deviation (Tci. This of course implies that 
even in the classical limit, m is not itself a Gaussian 
variable. 

The above behaviour is easily understood. Writing the 
complex cavity amplitude a{t) in terms of its mean value 
yH^g-icjrft g^j-^j ^ thermally fluctuating part 5a{t), we 
have: 



m{t) = / dt' 

'0 

dt' 



Vndr + e*"''*''5a(t') 







5X{t) + i5Y{t) 



(18) 



In the second line, we have written the fluctuation 5a{t) 
in terms of real- valued quadratures 5X{t) and 5Y{t). In 
the large fidr limit, only the low- frequency part of the 
"intensity-quadrature" noise X will be enhanced by fidr. 
To a good approximation, we may thus drop the Y con- 
tribution, and replace X by its low frequency part. We 
thus have: 



m(t) 



1 



t X 2 

dt'5X{t') 



{VW^t + 6X{t)) 



(19) 



One can easily confirm that the time-averaged intensity- 
quadrature noise 6X(t) defined above is a Gaussian ran- 
dom variable with zero mean and variance 2nthx/l'i 



Eq. (19) is thus in agreement with Eq. (16). It is worth 
noting that a full classical calculation of P{m) (includ- 
ing purely thermal effects) yields an answer in complete 
agreement with the classical limit of the moment gener- 
ating functions given in Eqs. (14a) and (14b). 



IV. 



ENERGY FLUCTUATION STATISTICS IN 
THE QUANTUM REGIME 

A. Basic results 



It is tempting to make a simple extrapolation of the 
classical energy statistics to the quantum regime. Again, 
in the large n^r limit it is the amplification of the ther- 



mal intensity-quadrature fluctuations SX{t) (cf. Eq. ( 19 )) 
which determine P{m)\ one might expect that the only 



difference in the quantum case is that these quadrature 
fluctuations are driven by both thermal and zero-point 
force noise. We would thus expect the distribution to 
again be given by Eq. (16), with the simple modifica- 
tion that the variance ad in Eq.(17) should be increased 



to include zero-point fluctuations via the substitution 
"th nth + 1/2. 

However, as already mentioned in the introduction, 
this is not the case. Instead, the full quantum moment 
generating function in Eq. ( 14a ) is related to the classical 
one (cf. Eq. ([T5|) by the simple substitution: 



Adr[fc;»^th] = Adr,ci [fc; ?^th -5> 



«cff 



ncs[k] 



nth 



47' 



(20a) 
(20b) 



Thus, one shifts nth both by the constant 1/2 (reflecting 
the inclusion of zero point force noise), as well as by an 
imaginary, fc-dependent term. 

The fc-dependence of rieff implies non-trivial quantum 
corrections to the third cumulant and higher involving 
a kind of feedback, whereby higher-order cumulants de- 
pend on the form of lower-order cumulants. In the case 
of the third cumulant ((to"^)), one finds: 



dn 



th 



47' 



(21) 



where we use ((to^))c1' to denote the naive expectation 
for the jth cumulant: the jth classical cumulant obtained 
from Eq. (15), with the substitution fith — nth + 1/2- 

Heuristically, this feedback of lower moments into 
higher moments is analogous to the situation in so-called 
"cascaded" Langevin approaches [El [14]. One could 
heuristically obtain the feedback term in Eq. (21) us- 



ing such an approach, starting with the assumption that 
the effective thermal occupation nth in the classical dis- 
tribution fluctuates in a way that is driven by (and hence 
correlated with) Sm. Assuming that the fluctuations of 
(5nth are slow compared to those of n (allowing a two-step 
averaging procedure), one obtains: 



dn 



th 



{Snthit)-Sm,{t)). 

(22) 



This recovers Eq. (21 ) if we take ((5nth(i)-i5'Ti(t)) = — 1/47 

a heuristic fashion, cascaded 
been used successfully to un- 
of current fluctuations in elec- 
Here, we stress that this pic- 
a fully quantum calculation, 
form of the third cumulant. 



While usually derived in 
Langevin approaches have 
derstand higher cumulants 
tronic conductors [131 [H] . 
ture emerges directly from 
Turning to the explicit 
evaluating Eq. (21) yields: 



// 3\\ n^j-t 2 

7 



24 (1 + 2nth) 



(23) 



The first term here is just the expected classical answer; 
it is always positive. The second term here is the non- 
trivial "feedback" quantum correction; as it involves the 



-4 -2 2 4 6 8 10 
(m-m) I ypjn^ 

FIG. 1: The distribution P{m) of integrated energy fluctua- 
tions m of a driven resonator, evaluated at a time t — IO/7 a 
drive detuning of 5 — IO7, and a driving force strength which 
yields an average number of cavity quanta fidr = 5. The var- 
ious curves correspond to different resonator temperatures: 
nth ~ 3 (black wide-dashed), nth ~ 0.25 (red, small-dashed) 
and fith ~ (solid blue). As discussed in the text, for large 
drive detunings and low temperatures, the distribution P(m) 
fails to be positive definite. 



second moment of the classical distribution, it is lower- 
order in the susceptibility x than the first term. As a 
result, this correction can make the skewness negative for 
a sufficiently small x, something that is impossible clas- 
sically. Further, in the limit of a strongly-detuned drive 
(i.e. X ^ 1), the quantum skewness has a much larger 
magnitude (by a factor 1 /x) than the corresponding clas- 
sical answer. The quantum "feedback" corrections simi- 
larly enhance all higher moments over the corresponding 
classical answer in the large detuning limit. 



B. Negative probabilities at large drive detuning 

The enhanced role of the non-trivial quantum correc- 
tions (arising from the k dependence of ricft[A:]) in the 
large detuning limit \6\ 3> 7 is even more apparent if 
one looks at the form of the full distribution P{m). One 
finds that at sufficiently low temperature and large de- 
tuning, these quantum corrections lead to P{m) becom- 
ing non-positive definite (see Fig. [I]) . This can be demon- 
strated analytically by just using the first four cumulants 
of P{m). The solution of the Hamburger moment prob- 
lem [23] is a necessary and sufficient set of conditions for a 
set of moments to correspond to a positive-definite proba- 
bility distribution. Letting Cj to denote the jth cumulant 

scaled by the variance (e.g. C3 = {{'m?)) l[{{'m'^))Y^'^), 
the lowest-order Hamburger positivity constraint involv- 
ing the third moment is: 



0.0 



0.5 



1.0 



1.5 



FIG. 2: The integrated total negativity of the distribution 
P(m) as a function of temperature nth, for a time t = IO/7 
and a driving force amplitude which yields ndr ~ 5. From 
left to right, the three curves correspond to drive detunings 
of S — 37, 5 — 5^ and 5 = 77. One sees that increasing the 
magnitude of the detuning causes the negativity to persist to 
higher temperatures. 



In the long time limit and for large detunings, C3 cx 
while C4 is independent of x- The above constraint is 
thus violated by P(m) for a sufficiently large detuning; 
in the large-t limit, the condition for violation becomes: 



(25) 



Heuristically, the quantum "feedback" contribution to 



the third moment (second term in Eq. (21)) has made 



it too large for P{m) to be positive definite [5T] . 

In Fig. w e plot the distribution P{m) as obtained 



from Eq. (14a 



[Csf < C4 + 2. 



(24) 



for a large detuning 5 = IO7 and for 
various bath temperatures; as the temperature is low- 
ered, the distribution fails to be positive for values of 
m > (to) . Fig. p]) plots the integrated negativity of the 
distribution, N\P] — — J dmP{m)6{—P{m)) as a func- 
tion of temperature. One clearly sees that increasing the 
drive detuning causes the negativity to emerge at pro- 
gressively higher temperatures. 

In the large detuning limit of interest, one can show an- 
alytically (see Appendix B) that it is indeed the anoma- 
lously large magnitude of the third cumulant (second 
term in Eq. (23)) which causes the distribution to be- 
come non-positive definite. As shown in Appendix B, 
the large value of the third cumulant causes the distribu- 
tion to have the form of an Airy function convolved with 
a Gaussian; the oscillations of the Airy function cause 
the P{m) to drop below for values of m > (m). 



V. INTERPRETATION 

The interpretation of negative quasi-probabilities of 
the kind found here (negative "counting statistics") was 
first given by Nazarov and Kindermann [12j . We begin by 
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quickly summarizing their findings, and then extend their 
interpretation to argue that negative counting statistics 
correspond to the same sort of non-classical behaviour 
tested by Leggett-Garg inequalities; in particular, they 
imply that a "macrorealistic" classical picture of the fluc- 
tuations of n{t) is not possible. 



A. Infinite mass detector 

Nazarov and Kindermann considered an alternate ide- 
alized setup for measuring m, where the detector is an 
infinitely heavy mass, and the quantity to be measured 
(in our case n) is linearly coupled to the position x of the 
detector: 



-Axfi. 



(26) 



The detector-oscillator interaction is turned on for a time 
t. Classically, the interaction would simply shift the de- 
tector momentum an amount Ant — Am, while (due to 
its infinite mass) its position would be unchanged. Corre- 
spondingly, one would expect that in the quantum case, 
the final momentum distribution of the detector mass 
would be a convolution of the detector's initial momen- 
tum distribution and the desired probability distribution 
P{m). The only additional complication is that there 
may be a backaction effect: the distribution P{ni) may 
itself (via i?int) depend on the value of the detector mass 
position X. In our case the interaction Hamiltonian in 
Eq. ( 26 ) implies that different values of x correspond to 



different values of the resonator frequency (and hence 
drive detuning 5). As a result, different values of x will 
lead to different distributions P{m\5 = + Ax). In- 
cluding this effect, one would then expect the follow- 
ing relation between the detector mass Wigner function 
W{x,p]t) before and after the interaction: 

W{x,p;t) = I dmP{m;5 = 5o+ Ax)W{x,p- Am]Q). 



(27) 

This relation was rigorously derived in Ref. [12], with 
P{m, S) being the usual Keldysh-ordered distribution we 
have been considering. 

Consider the case where the detector mass is initially in 
a Gaussian state with zero means, a momentum variance 
CTp ~ Aaimp and a position variance ax = ha^x/A. We 
could then use the final momentum distribution of the 
mass to infer the distribution P{m) as: 



-Pmcas("^;^) = 



A I dxW[x,m/A;t] 
1 



' (m — m')^ 



exp 



dm' I dS'P[m'; S'] x 



■ exp 



24a 



The above equation provides us with an unambiguous 
way to interpret P{m). It tells us that the Keldysh or- 
dered P(m, 5) should be regarded as the "intrinsic" dis- 
tribution describing the fluctuations of m for a given 
fixed value of ^. In contrast, the measured distribu- 
tion -Pmoas('7i, (5) is this intrinsic distribution corrupted 
by the addition of measurement uncertainty. There is 
both an imprecision uncertainty CTimp in m coming from 
the momentum uncertainty in the detector mass initial 
state, and a backaction uncertainty ctba in ^ coming from 
the detector mass position uncertainty. The Heiscnbcrg 
uncertainty principle applied to the detector mass im- 
plies that (TinipCTBA > 1/2: one can never eliminate both 
these sources of measurement uncertainty. Nonethe- 
less, Eq. (28) gives us a way to deflne the underlying. 



measurement- noise free distribution Pirn). 

Turning to the issue of positivity, we must have that 
the measured distribution Pmcas('7i, 5) is positive definite, 
as it is just the flnal momentum distribution of the test 
mass. In the special case where P{m) is independent of 
X (i.e. a true QND measurement where there is no back- 
action), this constraint immediately implies that P{m) 
be positive definite. However, in the case relevant here, 
where backaction is important (i.e. different values of x 
affect the fluctuations of m), this is no longer required: 
P{m) can exhibit negativity in these cases without vio- 
lating the positivity of Pmoas("T-)- This is precisely what 
we we find in the driven number-fluctuation statistics at 
low temperature and large drive detuning. 

For our driven cavity, P{m, S) will only change appre- 
ciably when 5 is varied an amount ^ 7. To avoid a size- 
able backaction, we would thus want ctba 1- The 
Heisenberg uncertainty principle then implies (Timp 3> 
1/7, which sets a limit to the scale Am of any negative 
regions in P{m,S). These constraints are indeed obeyed 
by our results. 



B. Significance of negative probabilities 



With Eq. ( 28 1 in hand, we can now view the failure of 



(28) 



P{m) to be positive definite as a clear manifestation of 
non-classical behaviour in our driven resonator. Classi- 
cally, we would naturally think of the fluctuations of m 
in terms of random trajectories n{t) and a correspond- 
ing distribution function. The failure of the Keldysh or- 
dered P{m) to be positive means that even in the most 
highly idealized measurement setups, we cannot interpret 
the measured statistics in this way: they do not corre- 
spond to having added measurement noise to an under- 
lying classical stochastic process. This interpretation can 
only be salvaged if one is willing to accept that the in- 
trinsic distribution function describing the fluctuations is 
not positive deflnite. 

Not surprisingly, negativity in P[m) and the corre- 
sponding non-classical behaviour only emerges at suffi- 
ciently low temperatures. More subtle however is the 
dependence on drive detuning: negativity only occurs at 
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a sufficiently large- magnitude detuning \S\, and is always 
enhanced by increasing \d\. On a basic level, this is con- 
sistent with Eq. (28 1, which tells us that P{m) can only 
be negative if there is a backaction effect associated with 
measuring m. In our case, this backaction effect vanishes 
to leading order when (5 = [9 . Hence, it is reasonable 
that obtaining negativity requires a non-zero detuning 5. 

We stress that this non-classicality discussed here is 
very different than that associated with a non-positive 
definite Wigner function; here, we are not characteriz- 
ing the instantaneous state of a system, but rather the 
time-integrated fluctuations of an observable. There is 
rather a much more natural connection to the kind of 
non-classical temporal correlations that lead to violations 
of Leggett-Garg inequalities (LGIs) |16]. These inequali- 
ties constrain temporal correlations of a given observable 
0{t) in any classical theory which satisfies: 

1. Macrorealism: 0{t) has a definite value at all times. 



2. Noninvasive measurability: 0{t) can be measured 
without any backaction disturbance that would al- 
ter its subsequent evolution. 



We note that two recent experiments have reported vio- 
lation of an LGI [2511261. 



In some sense, the non-classicality associated with the 
negativity of P(m) is stronger than that associated with 
the violation of an LGI. Violating an LGI could simply be 
interpreted as indicating that the measurement is indeed 
invasive (i.e. there is backaction), without necessarily in- 
dicating a failure of macrorealism. In contrast, negativity 
of P{m) tells us more than simply "backaction exists". 
It tells us that a natural way of including backaction ef- 
fects classically (as additional measurement noise which 
smears an intrinsic probability distribution, c.f. Eq. (28)) 
is impossible. Note that in our system, backaction ef- 
fects remain important when measuring P(m) at non- 
zero detuning even in the more classical regime of high- 
temperatures; nonetheless, the distribution exhibits no 
negativity here. 

It also interesting to note that while standard LGIs 
involve two-time correlation functions, the non-classical 
behaviour found here only manifests itself when one con- 
siders higher-order correlation functions. Recall that the 
second moment given in Eq. ([7]) has a completely clas- 
sical form, where the third moment and higher have 
non-trivial quantum corrections stemming from the k- 
dependence of the effective thermal occupancy factor 
neff[A;]. Finally, the standard violation of an LGI involves 
a qubit undergoing Larmor oscillations [TB] . Similarly, in 
our driven resonator negative probabilities only emerge 
for large-magnitude drive detunings \5\ ^ 7, a regime 
where correlation functions of n{t) have a strong oscilla- 
tory behaviour. 



VI. MEASUREMENT ISSUES 

We end with a discussion of how one might experi- 
mentally detect evidence of the non-classical photon and 
phonon number fluctuations described in this paper. 



A. Qubit plus resonator measurement 

One approach would be to experimentally implement 
the model of Sec. II C[ where a qubit is coupled dis- 
persively to the photon number operator n of a driven 
cavity (c.f. Eq.(|8])). The evolution of the qubit phase 
(i.e. its off-diagonal density matrix element p^^) at var- 
ious values of the dispersive coupling k directly yields 
the cumulant generating function of the desired distri- 
bution P{m) (c.f. Eq. (10)). Such a measurement could 



be contemplated in a cavity QED or circuit QED sys- 
tems, where a two-level system (atom or superconduct- 
ing qubit) is coupled to an electromagnetic cavity. The 
required phase information could be extracted by using 
a Ramsey-interference technique, similar to the seminal 
experiments of Refs. [3l |4j. Unlike those experiments, 
the focus here is very different: the goal is learn about 
the way the cavity photon number fluctuates over a time 
t ^ 1/7, as opposed to probe its instantaneous value. 

In order to detect evidence of the non-classical photon 
number fluctuations discussed here, it would be sufficient 
to see that the third moment {{m?)) is negative. Given 
the dispersive qubit-cavity interaction, the odd moments 
of P{m) will contribute in the long-time limit t':$> f/7 to 
the ac-Stark shift Ailqb of the qubit frequency. One has: 



qb 



lim I k 

t— >-C30 



knr\ 



m) 
~t 
^3 



fc3 ((m3)) 
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0(fc5) 



Oik^'xl 



(29) 



In the second line, we have assumed that the driven cav- 
ity is in the interesting regime of zero temperature and 
strong detuning (i.e. x ^ 1) where we expect strong 
quantum effects. We see that anomalous negative skew- 
ness of P{m) manifests itself in the sign of the k^ term 
in the ac-Stark shift of the qubit frequency. 

To resolve this non-linear contribution to the ac-Stark 
shift, one requires sufficient qubit coherence. As we must 
allow the qubit phase to evolve long enough both to be 
in the long-time limit of the photon-number fluctuations, 
and to resolve the k^ stark shift, we need that the total 
dephasing rate (including the contribution from Ti pro- 
cesses) satisfy: 



(30) 



< min ( 7, ^^^drX 



One unavoidable contribution to F^ will come from the 
dispersive coupling and the even moments of P(m); it 
is easy to see that this contribution satisfies Eq. (30) as 
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long as one has a strong dispersive coupling ^ 7. This 
regime has been achieved in several recent circuit QED 
experiments [SJ One also needs the intrinsic, non- 
cavity dephasing of the qubit to be sufficiently small. 
Given recent advances in extending the coherence of su- 
perconducting qubits [27 , this also would appear to be 
feasible. 

An alternate approach would be to measure the or- 
der term in the ac-Stark shift via simple spectroscopy, 
where one directly drives the qubit and measures its state 
as a function of this drive frequency. In order to avoid 
complications arising from the spectroscopy drive modi- 
fying the cavity state (and hence P{m)), one would need 
to use, e.g., a second cavity for the spectroscopy [5]. 



form: 



(("^'»BA 



1 



72 (1-1-4(52/72)4 



<32) 



By choosing a sufficiently small measurement strength 
(i.e. Sff) and large enough detuning, one can still have 
the total third moment be negative. In the large detun- 
ing limit, the backaction-induced skewness scales as \/5^, 
whereas the intrinsic, negative skewness scales as l/S"^. 
Note also that the backaction contribution to the second 
moment {{m?)) for this generic linear-detector setup was 
discussed in Ref. jjj; there, one finds that the minimum 
possible total added noise is achieved for Sff l/ndr- 



VII. CONCLUSIONS 



B. Measurement with a generic linear-response 
detector 



Perhaps a more general way to measure the fluctua- 
tions of m would be to weakly couple the photon (or 
phonon) number operator fi we wish to measure to the 
input port of a generic linear detector fT^ , as discussed 
in Ref. [9] . We would thus have a coupling of the form: 



int.lii 



hh ■ F 



(31) 



where -F is a detector operator and a generalized force. 
We would then monitor some other detector observ- 
able, say /, whose value depends linearly on fi. The 
dispersively-coupled optomechanical setup for detecting 
phonon shot noise analyzed in Ref. [Hj falls into this gen- 
eral description. In that case, F is the photon number 
operator of the optical cavity used to detect mechanical 
quanta, and / is the homodyne current. 

In this sort of generic setup, the statistics of the detec- 
tor output / can be used to to extract the statistics of m. 
Of course (similar to the idealized detector of Sec. VA|, 
this correspondence will be hindered by the presence of 
both measurement imprecision noise (i.e. the intrinsic 
fluctuations in /) as well as backaction noise (i.e. the 
effective fluctuations in detuning resulting from fluctua- 
tions in F). The simplest evidence for non-classical fluc- 
tuations of m comes from the anomalous sign of the third 
moment; we thus need to ask whether measurement im- 
precision and backaction would obscure the intrinsic neg- 
ativity of the skewness. 

To that end, we first note that measurement impreci- 
sion here can be treated as an additive Gaussian noise 
process, and hence will not change the value of the third 
moment. As for the backaction fluctuations, they are 
equivalent to having phase fluctuations on the mechan- 
ical drive. Treating these backaction phase fluctuations 
along the same lines as Ref. |28], we flnd that they yield 
an additional additive contribution to {{rrv^)) which is al- 
ways positive, and which in the large ridr limit takes the 



In this paper, we have shown that the full-counting 
statistics of energy fluctuations in a driven quantum res- 
onator can become negative for sufficiently low temper- 
ature and large drive detuning. This negativity arises 
from the same kind of quantum correction that leads to 
a negative third moment j9j, something that is impos- 
sible classically. We have argued that the failure of the 
quasi-probability distribution describing P{m) to be pos- 
itive definite is similar to having violated a Leggett-Garg 
inequality, and implies that a "macrorealistic" , classical 
picture for the fluctuations of phonon/photon number is 
not possible. We have thus shown in a relatively simple 
setting that higher moments of such counting statistics 
can be used to detect non-classical behaviour. It would 
be extremely interesting to investigate whether similar 
effects manifest themselves in other system. 
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Appendix A: Keldysh ordering from homodyne 
measurement theory 

Consider the measurement setup of recent optome- 
chanical experiments [71 [S], where the number opera- 
tor h{t) of a mechanical mode is coupled dispersively 
(strength g) to a driven optical mode which acts as a 
detector cavity. The coupling takes the form: 



Hi, 



gha) a. 



(Al) 



where a is the annihilation operator for the measurement 
cavity. By virtue of this interaction, the frequency of the 
detector cavity will depend on the value of n{t). One can 
thus measure the time variation of n{t) by detecting the 
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resulting variation in the detector cavity frequency via 
homodyne detection of its output field (see, e.g., [I9 1 I29 ] ). 
This involves first mixing the detector cavity output field 
with a large, classical reference beam. To leading order 
in g, the output field b from the mixer will have the form: 



b{t) = P + Bh(t). 



(A2) 



where [3 parameterizes the large magnitude of the clas- 
sical reference beam used, B on g and we have omitted 
vacuum noise terms responsible for shot noise (i.e. the 
imprecision noise in this measurement scheme) . Without 
loss of generality, wc take both f3 and B to be real. The 
intensity I — b^b oi the mixer output is then measured 
using a photodetector. Assuming that the constant term 
intensity can be subtracted from / (by, e.g., using bal- 
anced homodyne detection), the output of the detector 
to leading order in /? is: 



51 - (3^ ^ 2Bl3n. 



(A3) 



We see that in the large /3 limit, the output is just linearly 
proportional to fi. 

Given this simple linear correspondence, one can di- 
rectly infer the value of the jth moment of m from the 
measured jth-order intensity correlation function, i.e. 



meas — lim 

/3— )-oo 



1 

(2/3B)J 



t 3 



Yldt\{f si{t[)...i{t'^) 



1=1 



(A4) 



On the RHS of this expression, the symbol T denotes that 
the measured intensity correlation functions correspond 
to expectation values which are both normal-ordered and 
time-ordered with respect to the b{t) and &^(t) operators; 
this ordering prescription is a direct consequence of mea- 
suring intensity via photodetection [2D1[3U]. For example, 
for the third moment we have: 

f{I{h)I{h)I{h)) = {b^ta)bHh)b\tSb{tc)b{h)b{ta)), 



where ta < h 
ti,t2 and ts. 



(A5) 

< tc denotes the time-ordered listing of 



Eq. ( A4 ) and Eq. ( A2 ) completely determine the corre- 
spondence between the measured moments (?Ti-')moas and 
appropriately ordered expectation values of products of 
m{t). Though tedious, one can now explicitly confirm 
that for each moment (m^)nicas, the resulting ordering of 
m{t) operators is exactly the Keldysh ordering defined by 
Eq. ( |lO| ). For the third moment, this was done in Ref. [S]. 

A more compact way to see that one obtains the 
Keldysh ordering at each order is to use a functional- 
integral formulation of the Keldysh technique; a peda- 
gogical introduction to this approach is given in Refs. [U] 
[25]. In this formulation, each bosonic operators is 
replaced by two time-dependent fields, e.g. b{t) — 
ba-{t),n{t) — >■ na-{t), where the index a — +{—) denotes 



the forward (backwards) Keldysh contour. Within this 
approach, different operator orderings correspond to dif- 
ferent combinations of + and — fields. A special role 
is played by the so-called "classical" field, which is the 
average of -|- and — fields, i.e. 



nci(i) 



(A6) 



At the level of a saddle-point approximation, the dynam- 
ics of the classical field correspond to an effective classical 
Langevin equation. Correlation functions of this classi- 
cal field are obtained in the usual way using the Keldysh 
action S describing the system: 



{nci{ti)...nci{tj)) = J Y[ 'D(f>j„nci{ti)...nci{tj)exp[iS]. 



(A7) 



The (j>ja{t) here denote the various fields that describe 
the system; the action is a function of these fields. By 
construction, the jth-order correlation functions defined 
above is identical to the jth order, Keldysh-ordered op- 
erator expectation value defined by Eq. (10) [21]. 



Turning to our homodyne measurement, we first note 
that the normal-ordered, time-ordered intensity correla- 
tion functions that are measured via photodetection can 
be obtained by adding an auxiliary source term to the 
Keldysh action S of the form: 



k I dt' [b*_{t)b+{t) ~ 13^ 



(A8) 



Given the correspondence between the Keldysh ± fields 
and operator orderings ^22j, one finds that derivatives of 
the full Keldysh partition function (action S + Ssrc) with 
respect to fc at fc = generate the desired normal and 
time-ordered correlation functions in the usual way. 

Next, for homodyne detection, we can make the re- 
placement: 



b^{t) ~ bl{t)~/3 + Bn,{t), 



which results in: 



2kf5B / dt'riciit), 
lo 



(A9) 



(AlO) 



where we have only retained the leading-order- in-/? term 
in the action. We thus see that the source field k cou- 
ples to the classical field nci{t); it thus follows that the 
jth-order intensity correlation functions (as generated by 
differentiation of the Keldysh partition function with re- 
spect to k) will be directly proportional to jth order 
Keldysh-ordered correlation functions of m. 



Appendix B: P{m) in the large time, large detuning 

limit 



We first shift and rescale P{m) in the full quantum 
case so that it has zero mean and unit variance. Setting 
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nth = 0, the CGF Adi [A:] of the transformed distribution 
takes the form: 



fc2 

^ » 1 + i 



1 



ndrt7 



(Bl) 



4ndri7 



Consider now the strong-detuning, long-time hmit, 
such that fidilt od but x^dijt is finite. In this limit 



1 



(B2) 



In the quantum case, both the second and third moments 
are non-vanishing in this long-time, strong-detuning 
limit. In contrast, in the same limit the classical dis- 
tribution would be completely Gaussian. Thus, the third 
moment term in Eq. (B2 1 is entirely due to the effective k 



dependence of the thermal factor fieff in the quantum dis- 
tribution. Further, note that increasing the drive detun- 
ing (and hence reducing x) enhances the non-Gaussian 
nature of the distribution described by Eq. ( B2 ) ; this is 



the opposite of what happens classically, where a large 
detuning suppresses fluctuations and non-Gaussian ef- 
fects, as the magnitude of thermal fluctuations at the 
drive frequency are reduced. 



Fourier transforming the approximate CGF in Eq. (B2 ) 
reveals that in the large detuning limit, the distribu- 



tion P{m) is a convolution of a Gaussian and an Airy 
function. This can be explicitly evaluated. Defining 

m = (to — fid^t) I ((to^))), we have: 



P{m,t) 



where 



1 



■ cxp 



-1 
2A3 



1 
6A3 



Ai 



TO 



{{^')) 3 

2((to2))3/2 ^ 8V^I^' 



1 

(B3) 



(B4) 



It is the oscillation of the Airy function factor in Eq. ( B3 ) 



above which gives rise to the negative probabilities at 
TO > 0. The exponential prefactor ensures that the re- 
sulting negativity is exponentially suppressed in the long 
time limit when A <C 1. However, for intermediate times 
(still much longer than I/7), one has A ^ 1, and the neg- 
ativity can be appreciable. Note that the most prominent 
domain of negativity in this large-A limit has an extent in 
TO ~ A; in terms of m/t, this corresponds to a range < 1. 
Thus, while P{m) exhibits negativity even in a seem- 
ingly classical regime where n^-c^t ^ 1, it only occurs on 
a scale which corresponds to less than one quantum in 
the resonator. This is consistent with the discussion of 
negative quasi-probabilities given in Sec. |V A| 
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